#import libraries
library(psych)
library(dplyr)
library(effectsize)
library(pscl)
library(pwr)
library(standardize)
library(Hmisc)
library(ltm)


#load data
data <- read.csv("Study 2.csv")
View(data)


#impute missing value of experience
##two participants input wrong information about their experience
##Impute their values with median as the distribution of experience is not normal
hist(data$exp_length)
data$exp_length[is.na(data$exp_length)] <- median(data$exp_length, na.rm = TRUE)



#recode categorical vars with values
data$gender_coded <- recode_factor(data$gender, 
                               '1' = 'male', 
                               '2' = 'female', 
                               '3' = 'non_binary_third_gender', 
                               '4' = 'other')

data$bad_response_coded <- recode_factor(data$bad_response, 
                                         '0' = 'passed', 
                                         '1' = 'failed')


#get the sample information
##gender
gender_count<-table(data$gender_coded) #71 males vs. 65 females vs. 2 binary/third-gender vs. 1 not to say
gender_count

##age
summary(data$age) #37.12
sd(data$age) #9.55

##bad respondents
bad_response_count<-table(data$bad_response_coded) #108 passed vs. 31 failed
bad_response_count


#calculate inequality as the focal IV
data <- data %>%
  mutate(inequality_mean = (inequality_1 + inequality_2 + inequality_3)/3)

#reliability test of inequality
inequality_alpha <- data[, c('inequality_1', 'inequality_2', 'inequality_3')]
alpha(inequality_alpha) #alpha = .77

#descriptive information of continuous focal IVs and DV
focal_cont_vars <- data[, c('willingness', 
                        'inequality_mean', 'inequality_1', 'inequality_2', 'inequality_3',
                        'exp_length','income')]

for (col_name in colnames(focal_cont_vars)) {
  col_data <- focal_cont_vars[[col_name]]
  col_mean <- mean(col_data, na.rm = TRUE)  
  col_sd <- sd(col_data, na.rm = TRUE)
  

  cat("Variable:", col_name, "\n")
  cat("Mean (M):", col_mean, "\n")
  cat("Standard Deviation (SD):", col_sd, "\n\n")
}
###Variable: willingness  Mean (M): 5.669065  Standard Deviation (SD): 1.293173 
###Variable: inequality_mean Mean (M): 5.067146  Standard Deviation (SD): 1.074653 


#correlation matrix of continuous focal IVs and DV
cor_matrix <- rcorr(as.matrix(focal_cont_vars))
cor_matrix$r #r = -.186
cor_matrix$P #p = .029

cor_test <- cor.test(data$inequality_mean, data$willingness,use = "complete.obs") #t(137)  = -2.21
cor_test

#regression
model_focal <- lm(willingness ~ inequality_mean, data = data)
summary(model_focal)
standardize_parameters(model_focal)

model_cov <- lm(willingness ~ inequality_mean + exp_length + income + gender + race_White + host_exp, data = data)
summary(model_cov)
standardize_parameters(model_cov)


#repeat the analysis for participants passed the check
selected_data <- data %>% filter(bad_response == 0)
View(selected_data)

focal_cont_vars_selected <- selected_data[, c('willingness', 
                            'inequality_mean', 'inequality_1', 'inequality_2', 'inequality_3',
                            'exp_length','income')]

inequality_alpha_selected <- selected_data[, c('inequality_1', 'inequality_2', 'inequality_3')]
cronbach.alpha(inequality_alpha_selected) #alpha = .79

for (col_name in colnames(focal_cont_vars_selected)) {
  col_data <- focal_cont_vars_selected[[col_name]]
  col_mean <- mean(col_data, na.rm = TRUE)  
  col_sd <- sd(col_data, na.rm = TRUE)
  
  
  cat("Variable:", col_name, "\n")
  cat("Mean (M):", col_mean, "\n")
  cat("Standard Deviation (SD):", col_sd, "\n\n")
}
###Variable: willingness Mean (M): 5.685185 Standard Deviation (SD): 1.357899 
###Variable: inequality_mean  Mean (M): 5.055556  Standard Deviation (SD): 1.082764 

cor_matrix_selected <- rcorr(as.matrix(focal_cont_vars_selected))
cor_matrix_selected$r #r = -.217
cor_matrix_selected$P #p = 0.024 

cor_test2 <- cor.test(selected_data$inequality_mean, selected_data$willingness,use = "complete.obs") #t(106)  = -2.29
cor_test2

model_focal_selected <- lm(willingness ~ inequality_mean, data = selected_data)
summary(model_focal_selected)
standardize_parameters(model_focal_selected)

model_cov_selected <- lm(willingness ~ inequality_mean + exp_length + income + gender + race_White + host_exp, data = selected_data)
summary(model_cov_selected)
standardize_parameters(model_cov_selected)